
if YesGPU==1
    PROBS=gpuArray(PROBS);
    Z=gpuArray(Z);
end

param=beta;
%param=horzcat(gamma, beta(1,:));
%for j=2:NCLASS
%    param = horzcat(param, beta(j,:));
%end

disp('Start estimation');
disp('The negative of the log-likelihood is minimized,');
disp('which is the same as maximizing the log-likelihood.');
tic;
% was 'iter' after 'Display'
options=optimset('LargeScale','off','Display','off','GradObj','on',...
    'MaxFunEvals',10000,'MaxIter',MAXITERS,'TolX',PARAMTOL,'TolFun',LLTOL,'DerivativeCheck','off');
if WantHessian==1;
    [paramhat,fval,exitflag,output,grad,hessian]=fminunc(@flexll,param,options);
else
    [paramhat,fval,exitflag]=fminunc(@flexll,param,options);
end
finalLL=-flexll(paramhat);
disp(' ');
disp(['Estimation took ' num2str(toc./60) ' minutes.']);
disp(' ');
disp('Calculating summary statistics for random utility parameters.');

if YesGPU==1;
    PROBS=gather(PROBS);
    Z=gather(Z);
end;

[MeanEst,StdEst,CovMatEst,FreqEst,MidEst]=stats(paramhat,NBins);
disp('Estimated coefficients of Z variables are held in paramhat.');
if WantHessian==1
    disp('Hessian at convergence is held in hessian.');
end;
disp('Means, Standard Deviations, and Covariance Matrix of Utility Parameters');
disp('are held as MeanEst, StdEst, and CovMatEst.');
disp('Share of density in each bin and midpoint for each bin are held');
disp('in FreqEst and MidEst; each is size NV x NBins');
disp('Means and StdDevs of Random Utility Paramaters');

disp(' ');
if CrossCorr==1
    disp('Correlation Matrix for WTPs');
    if discountfun~="QuasiHyper"
        disp(corrcov(CovMatEst(3:end,3:end)));
    else 
        disp(corrcov(CovMatEst(4:end,4:end)));
    end
end
disp(' ');
disp('                Mean     Std Dev');
disp('              -------------------');
for r=1:length(NAMES);
    fprintf('%-10s %10.4f %10.4f\n', NAMES{r,1}, [MeanEst(r,1) StdEst(r,1)]);
end
disp(' ');
disp(' ');
disp('Final Log-Likelihood');
disp(finalLL);
disp(' ');
disp('Final AIC');
finalAIC = 2*NZ -2*finalLL;
disp(finalAIC);
disp(' ');
disp('Final BIC');
finalBIC = NZ*log(NP*5) - 2*finalLL;
disp(finalBIC);
disp(' ');

%disp('To obtain histogram for utility parameter k, type: bar(MidEst(k,:),FreqEst(k,:))');


